Quality metrics per cell

colData(sce)$total_counts[1:5] #Total count for each cell

## [1] 29060 29304 38929 40557 27625

colData(sce)$total_features_by_counts[1:5] #Number of genes with non zero counts per cell

## [1] 4898 4750 6090 5887 4753

colData(sce)$pct_counts_Spike[1:5] #Percentage of spike-in counts per cell

## [1] 23.07639 21.95946 16.27322 17.33856 21.46968

Quality metrics per gene

rowData(sce)$mean_counts[1:5] #Average count per gene

## [1] 0.27886855 0.42362729 1.04292845 0.04459235 0.37936772

rowData(sce)$n_cells_by_counts[1:5] #Number of cells with non zero counts per gene

## [1]  472  573 1167   77  669

rowData(sce)$pct_dropout_by_counts[1:5] #Percentage of dropouts per gene

## [1] 84.29285 80.93178 61.16473 97.43760 77.73710

Plot histograms based on quality metrics

One can plot histograms based on quality metrics to determine possible thresholds for cell/gene filterings.

hist(sce$total_counts, breaks=50,  xlab="Library size", ylab="Number of cells")

hist(sce$total_features_by_counts, xlab="Number of expressed genes", breaks=50,  ylab="Number of cells")

hist(sce$pct_counts_Spike, xlab="ERCC proportion (%)", ylab="Number of cells", breaks=50)

Automatic filtering of outlier cells based on quality metrics

Alternatively, isOutlier function provides an automatic detection of low quality cells based on quality metrics. In this example, we are filtering cells based on library size, number of expressed genes per cell and total count over all spike-in transcripts in each cell. Excluded cells are those with the total number of expressed genes and the total sum of counts more than 3 median absolute deviations below the median across the genes.

libsize.drop <- isOutlier(sce$total_counts, nmads=3, type="lower", log=TRUE)
feature.drop <- isOutlier(sce$total_features_by_counts, nmads=3, type="lower", log=TRUE)
spike.drop <- isOutlier(sce$pct_counts_Spike, nmads=3, type="higher")
sce <- sce[,!(libsize.drop | feature.drop | spike.drop)]
data.frame(ByLibSize=sum(libsize.drop), ByFeature=sum(feature.drop), BySpike=sum(spike.drop), Remaining=ncol(sce))

##   ByLibSize ByFeature BySpike Remaining
## 1         8         3       8      2989

Explore highly expressed genes

We can also explore highly and overexpressed genes. The overexpressed genes, are those genes which expression is higher than others by several magnitudes. Overexpressed genes can bias further procedures such as clustering, therefore, they should be removed before performing downstream analysis. Note that highly expressed genes mostly include spike-ins and mitochondrial genes. One of the most commonly overexpressed endogenous gene is “Malat1”.

plotHighestExprs(sce, exprs_values = "counts")

Filter overexpressed features

In this step we will filter spike-in and mitochondrial genes as well as “Malat1”.

sce=sce[-which(grepl("^ERCC-", rownames(sce))),]
sce=sce[-which(grepl("^mt-", rownames(sce))),]
sce <- sce[-which(rownames(sce)=="Malat1"),]
dim(sce)

## [1] 19804  2989

Filter lowly expressed genes based on the average count per gene

Now, we should filter out lowly expressed genes as they do not provide any insights into the underlying biology. In the filtering we removed lowly expressed genes that are genes with average expression count (adjusted by library size) equal to 0.

rowData(sce)$ave.counts <- calcAverage(sce, exprs_values = "counts", use_size_factors=FALSE)
to.keep <- rowData(sce)$ave.counts > 0
sce <- sce[to.keep,]
summary(to.keep)

##    Mode   FALSE    TRUE 
## logical       2   19802

dim(sce)

## [1] 19802  2989

Perform standard normalization by computing counts per million (CPM)

After quality control and filtering, one should normalize the dataset to remove potential technical bias. One strategy is to use CPM (Count Per Million) normalization which is a correction to remove the noise related to sequencing depth. CPM divides each count by its total sum (across all the genes) and multiplies by one million. In this way, each cell has the same total sum of the counts.

cpm(sce) <- calculateCPM(sce, use_size_factors=FALSE)
cpm(sce)[1:4,1:4]

##                V3       V4      V5        V6
## Tspan12    0.0000  0.00000   0.000  94.75081
## Tshz1    144.5226 47.89501   0.000  63.16720
## Fnbp1l   144.5226 47.89501 197.336 126.33441
## Adamts15   0.0000  0.00000   0.000   0.00000

Perform normalization using decnvolution method

Another strategy, more suitable for scRNAseq data, is deconvolution method implemented in scran package. The deconvolution method normalizes data by cells-pooled size factors that account for dropout biases. More details about this normalization technique can be found in study https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0947-7.

set.seed(100)
clusters <- quickCluster(sce, min.size=200, min.mean=0.1, method="igraph")
sce <- computeSumFactors(sce, cluster=clusters, min.mean=0.1)
sce <- scater::normalize(sce, return_log = FALSE)
normcounts(sce)[1:4,1:4]

##                V3        V4       V5       V6
## Tspan12  0.000000 0.0000000 0.000000 1.278561
## Tshz1    1.894734 0.6425685 0.000000 0.852374
## Fnbp1l   1.894734 0.6425685 2.582446 1.704748
## Adamts15 0.000000 0.0000000 0.000000 0.000000

Feature selection to retrieve most variable genes

Feature selection step is aimed at preserving biologically relevant information in the dataset and improving computational efficiency of downstream analyses. In most of the cases, it seeks for highly variable genes based on mean/variance relationship. M3Drop is one of the packages for automatic detection of most variable genes. It first calculates the mean and square coefficient of variation and fits the quadratic curve between the two variables. Then a chi-square test is used to find high variable genes which are significantly above the curve. We found CPM adjusted counts more appriopriate to use as input for this method rather than scran normalized values.

hvg_genes <- BrenneckeGetVariableGenes(cpm(sce), fdr = 0.01, minBiolDisp = 0.5)

length(hvg_genes)

## [1] 8183

sce_sub=sce[hvg_genes,]

Dimension reduction to project data into low-dimensional embedding

Dimension reduction can be used to visualize basic structure of the data or to reduce the number of features prior downstream analysis such as clustering. In contrast to the feature selection which extracts the most informative features, dimension reduction techniques projects the data into a new low-dimensional space that preserves the structure of the data. To reduce dataset dimensionality one can use PCA, tSNE or UMAP techniques implemented in scater package. Note that projections should be applied to normalized and log-transformed count matrices. PCA is a deterministic approach while tSNE and UMAP are stochastic - for reproducibility of tSNE and UMAP projections one has to set the seed for generating random variables.

assay(sce, "logcounts") <- log2(normcounts(sce)+1) 
plotPCA(sce, colour_by="level1class", run_args=c(exprs_values = "logcounts"))  

## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.

set.seed(100)
plotTSNE(sce, colour_by="level1class", run_args=c(exprs_values = "logcounts"))

## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.

plotUMAP(sce, colour_by="level1class", run_args=c(exprs_values = "logcounts"))

## Warning: 'add_ticks' is deprecated.
## Use '+ geom_rug(...)' instead.

Exercise

How the projections would change if you would plot raw or CPM counts followed by log-transformation?

Evaluate clustering output based on annotation

If the cell annotation is available, one can measure the effectiveness of clustering output using Adjusted Rand Index (ARI). ARI index measures similarity between the partition obtained from clustering and the partition obtained from dataset annotation. The values of the ARI range can be negative if the agreement of the partitions is worse then the agreement expected by chance, or between 0 and 1 for clustering better then chance. ARI close to 1 indicate high accuracy of the method in detecting annotated cell populations.

adjustedRandIndex(colData(sce)[paste0("sc3_",k_input, "_clusters")][[1]],sce$level1class)

## [1] 0.3084636

Exercise

How the clustering would change if you would use true number of cell populations (taken from annotation) in the inference of SC3?

Create Seurat object

seur <- CreateSeuratObject(raw.data = counts(sce), meta.data = data.frame(colData(sce)), min.cells = 0, min.genes = 0)
seur

## An object of class seurat in project SeuratProject 
##  19802 genes across 2989 samples.

Seurat also provides function to normalize data by scaling factor followed by log-transformation.

seur <- NormalizeData(seur, normalization.method = "LogNormalize")
seur@data[1:4,1:4]

## 4 x 4 sparse Matrix of class "dgCMatrix"
##                 V3        V4       V5        V6
## Tspan12  .         .         .        0.6665506
## Tshz1    0.8941375 0.3913325 .        0.4896053
## Fnbp1l   0.8941375 0.3913325 1.089693 0.8168434
## Adamts15 .         .         .        .

Extract highly variable genes for further analysis

For extracting most informative genes one can use FindVariableGenes function. Similarly to other feature selection strategies FindVariableGenes selects informative features beased on mean / dispersion relation.

seur <- FindVariableGenes(object = seur)

Scale the data

Scales and centers genes in the dataset.

seur<- Seurat::ScaleData(object = seur)

## Scaling data matrix

Reduce dataset dimension with PCA and clusterize cells

PCA or ICA (Independent Components Analysis), both with a user-specified or the internal number of reduced dimensions () (the internal value is equal to the full number of features). Then, on the data reduced space, K-nearest neighbor (KNN) graph is constructed (based on the Euclidean distances between any pair of cells) and weighted by the shared overlap in the local neighborhoods of the nodes. to cluster cells represented in such graph structure, uses modularity detection Louvain algorithm (herein called ) that estimates the number of clusters by optimizing the modularity function .

seur <- RunPCA(object = seur, pc.genes=seur@var.genes)
seur <- FindClusters(object = seur, reduction.type="pca")

table(seur@ident)

## 
##   0   1   2   3   4   5   6   7 
## 853 666 372 311 276 246 183  82

Compare clusterization with annotated cell types

Seurat::DimPlot(seur, group.by = "ident", reduction.use = "pca")

Seurat::DimPlot(seur, group.by = "level1class", reduction.use = "pca")

adjustedRandIndex(as.numeric(seur@ident), seur@meta.data$level1class)

## [1] 0.7844589

Create CellDataSet data object

CellDataSet data object is similar for SingleCellExperiment object. It stores information about the experiment (i.e. expression matrix and metadata) in “slots”. Note that when creating a new CellDataSet data should not be normalized - Monocle, when performing downstream analysis steps, normalizes the data internally.

rowData(sce)$gene_short_name <- rownames(sce)
pd <- new("AnnotatedDataFrame", data = data.frame(colData(sce))) #cell info
fd <- new("AnnotatedDataFrame", data = data.frame(rowData(sce))) #gene info
cds <- newCellDataSet(cellData = counts(sce), phenoData = pd, featureData = fd, expressionFamily=VGAM::negbinomial.size())  
cds

## CellDataSet (storageMode: environment)
## assayData: 19802 features, 2989 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: V3 V4 ... V3007 (2989 total)
##   varLabels: tissue group ... Size_Factor (58 total)
##   varMetadata: labelDescription
## featureData
##   featureNames: Tspan12 Tshz1 ... Gm21943 (19802 total)
##   fvarLabels: is_feature_control is_feature_control_Spike ...
##     gene_short_name (17 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:

Extract highly variable genes for further analysis

cds <- BiocGenerics::estimateSizeFactors(cds)
cds <- estimateDispersions(cds)

## Removing 129 outliers

disp <- dispersionTable(cds)
head(disp)

##    gene_id mean_expression dispersion_fit dispersion_empirical
## 1  Tspan12      0.27716042       5.633211            10.023476
## 2    Tshz1      0.41719267       3.973567             8.826994
## 3   Fnbp1l      0.83724181       2.325529             3.491744
## 4 Adamts15      0.04898606      28.664478            83.258508
## 5   Cldn12      0.30621583       5.164049             5.294964
## 6    Rxfp1      0.03631543      38.425363            12.890793

hvg <- subset(disp, mean_expression >= 0.5 & dispersion_empirical >= 0.1)
cds <- setOrderingFilter(cds, ordering_genes = hvg$gene_id)
plot_ordering_genes(cds)

## Warning: Transformation introduced infinite values in continuous y-axis

cds <- cds[hvg$gene_id,]

Reduce dimension of the data and cluster cells to detect cell populations

cds <- reduceDimension(cds, max_components = 2, reduction_method = 'tSNE', verbose = T)

## Warning in if (method == "PCA") {: the condition has length > 1 and only
## the first element will be used

cds <- clusterCells(cds, num_clusters = NULL)

## Distance cutoff calculated to 3.26211

Compare clusterization with annotated cell types

p1 <- monocle::plot_cell_clusters(cds, color_by = 'Cluster')
p1

p2 <- monocle::plot_cell_clusters(cds, color_by = 'level1class')
p2

adjustedRandIndex(cds$Cluster,cds$level1class)

## [1] 0.6961306

Package versions used in this analysis

sessionInfo()

## R version 3.5.1 (2018-07-02)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
## 
## Matrix products: default
## BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
## LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] SC3_1.10.1                  scran_1.10.2               
##  [3] mclust_5.4.3                Seurat_2.3.4               
##  [5] cowplot_0.9.4               monocle_2.99.2             
##  [7] L1Graph_0.1.1               lpSolveAPI_5.5.2.0-17      
##  [9] DDRTree_0.1.5               irlba_2.3.3                
## [11] igraph_1.2.4                Matrix_1.2-17              
## [13] M3Drop_1.8.1                numDeriv_2016.8-1          
## [15] scater_1.10.1               ggplot2_3.1.0              
## [17] SingleCellExperiment_1.4.1  SummarizedExperiment_1.12.0
## [19] DelayedArray_0.8.0          BiocParallel_1.16.6        
## [21] matrixStats_0.54.0          Biobase_2.42.0             
## [23] GenomicRanges_1.34.0        GenomeInfoDb_1.18.2        
## [25] IRanges_2.16.0              S4Vectors_0.20.1           
## [27] BiocGenerics_0.28.0        
## 
## loaded via a namespace (and not attached):
##   [1] prabclus_2.2-7           R.methodsS3_1.7.1       
##   [3] coda_0.19-2              pkgmaker_0.27           
##   [5] tidyr_0.8.3              acepack_1.4.1           
##   [7] bit64_0.9-7              knitr_1.22              
##   [9] R.utils_2.8.0            data.table_1.12.0       
##  [11] rpart_4.1-13             RCurl_1.95-4.12         
##  [13] doParallel_1.0.14        metap_1.1               
##  [15] snow_0.4-3               RANN_2.6.1              
##  [17] VGAM_1.1-1               combinat_0.0-8          
##  [19] proxy_0.4-23             future_1.12.0           
##  [21] bit_1.1-14               webshot_0.5.1           
##  [23] httpuv_1.5.0             assertthat_0.2.1        
##  [25] viridis_0.5.1            xfun_0.5                
##  [27] evaluate_0.13            promises_1.0.1          
##  [29] DEoptimR_1.0-8           caTools_1.17.1.2        
##  [31] DBI_1.0.0                htmlwidgets_1.3         
##  [33] sparsesvd_0.1-4          spdep_1.0-2             
##  [35] purrr_0.3.2              RSpectra_0.13-1         
##  [37] crosstalk_1.0.0          dplyr_0.8.0.1           
##  [39] backports_1.1.3          trimcluster_0.1-2.1     
##  [41] gbRd_0.4-11              deldir_0.1-16           
##  [43] Cairo_1.5-10             ROCR_1.0-7              
##  [45] withr_2.1.2              robustbase_0.93-4       
##  [47] checkmate_1.9.1          cluster_2.0.7-1         
##  [49] ape_5.3                  segmented_0.5-3.0       
##  [51] lazyeval_0.2.2           crayon_1.3.4            
##  [53] hdf5r_1.1.1              labeling_0.3            
##  [55] glmnet_2.0-16            edgeR_3.24.3            
##  [57] pkgconfig_2.0.2          slam_0.1-45             
##  [59] units_0.6-2              nlme_3.1-137            
##  [61] vipor_0.4.5              nnet_7.3-12             
##  [63] rlang_0.4.0              globals_0.12.4          
##  [65] diptest_0.75-7           miniUI_0.1.1.1          
##  [67] registry_0.5-1           doSNOW_1.0.16           
##  [69] ggrastr_0.1.7            lmtest_0.9-36           
##  [71] rngtools_1.3.1           Rhdf5lib_1.4.3          
##  [73] boot_1.3-20              zoo_1.8-5               
##  [75] base64enc_0.1-3          beeswarm_0.2.3          
##  [77] ggridges_0.5.1           pheatmap_1.0.12         
##  [79] png_0.1-7                viridisLite_0.3.0       
##  [81] bitops_1.0-6             R.oo_1.22.0             
##  [83] KernSmooth_2.23-15       DelayedMatrixStats_1.4.0
##  [85] rgl_0.100.19             doRNG_1.7.1             
##  [87] classInt_0.3-1           lars_1.2                
##  [89] stringr_1.4.0            manipulateWidget_0.10.0 
##  [91] scales_1.0.0             magrittr_1.5            
##  [93] plyr_1.8.4               ica_1.0-2               
##  [95] gplots_3.0.1.1           bibtex_0.4.2            
##  [97] gdata_2.18.0             zlibbioc_1.28.0         
##  [99] compiler_3.5.1           HSMMSingleCell_1.2.0    
## [101] lsei_1.2-0               bbmle_1.0.20            
## [103] RColorBrewer_1.1-2       rrcov_1.4-7             
## [105] fitdistrplus_1.0-14      dtw_1.20-1              
## [107] XVector_0.22.0           LearnBayes_2.15.1       
## [109] listenv_0.7.0            pbapply_1.4-0           
## [111] htmlTable_1.13.1         Formula_1.2-3           
## [113] MASS_7.3-51.3            tidyselect_0.2.5        
## [115] stringi_1.4.3            densityClust_0.3        
## [117] yaml_2.2.0               locfit_1.5-9.1          
## [119] latticeExtra_0.6-28      ggrepel_0.8.0           
## [121] pbmcapply_1.3.1          grid_3.5.1              
## [123] tools_3.5.1              rstudioapi_0.10         
## [125] foreach_1.4.4            foreign_0.8-71          
## [127] gridExtra_2.3            Rtsne_0.15              
## [129] digest_0.6.18            FNN_1.1.3               
## [131] shiny_1.2.0              qlcMatrix_0.9.7         
## [133] fpc_2.1-11.1             Rcpp_1.0.1              
## [135] SDMTools_1.1-221         later_0.8.0             
## [137] WriteXLS_4.1.0           httr_1.4.0              
## [139] sf_0.7-3                 npsurv_0.4-0            
## [141] kernlab_0.9-27           Rdpack_0.10-1           
## [143] colorspace_1.4-1         reticulate_1.11.1       
## [145] umap_0.2.0.0             splines_3.5.1           
## [147] statmod_1.4.30           expm_0.999-4            
## [149] sp_1.3-1                 flexmix_2.3-15          
## [151] plotly_4.8.0             spData_0.3.0            
## [153] xtable_1.8-3             jsonlite_1.6            
## [155] dynamicTreeCut_1.63-1    modeltools_0.2-22       
## [157] R6_2.4.0                 gmodels_2.18.1          
## [159] Hmisc_4.2-0              pillar_1.4.2            
## [161] htmltools_0.3.6          mime_0.6                
## [163] glue_1.3.1               BiocNeighbors_1.0.0     
## [165] class_7.3-15             codetools_0.2-16        
## [167] pcaPP_1.9-73             tsne_0.1-3              
## [169] mvtnorm_1.0-10           lattice_0.20-38         
## [171] tibble_2.1.1             mixtools_1.1.0          
## [173] ggbeeswarm_0.6.0         gtools_3.8.1            
## [175] survival_2.44-1.1        limma_3.38.3            
## [177] rmarkdown_1.12           docopt_0.6.1            
## [179] fastICA_1.2-1            munsell_0.5.0           
## [181] e1071_1.7-1              rhdf5_2.26.2            
## [183] GenomeInfoDbData_1.2.0   iterators_1.0.10        
## [185] HDF5Array_1.10.1         reshape2_1.4.3          
## [187] gtable_0.3.0